Supplementary Report S2. Report with ASV without OTU post-clustering. Terroir and farming practices drive arbuscular mycorrhizal fungal communities in French vineyard.

Code
knitr::opts_chunk$set(
  message = FALSE
)
Code
#if (!require("devtools", quietly = TRUE)) {
#  install.packages("devtools")
#}
#devtools::install_github("adrientaudiere/MiscMetabar")
Code
library(MiscMetabar)
library(targets)

library(metacoder)
library(DT)
library(knitr)
library(ggplot2)
library(patchwork)
library(vegan)
library(FactoMineR)
library(factoextra)
library(DESeq2)
library(indicspecies)
library(sf)
library(adespatial)
library(ade4)
library(adegraphics)
library(kableExtra)
library(formattable)
library(microbiomeMarker)
library(ComplexUpset)
Code
d <- tar_read(d_blast)

Metadata modification

Code
d@sam_data$rank <- gsub("Weeding", "Grassing", d@sam_data$rank)
d@sam_data$inter_rank <- gsub("Weeding", "Grassing", d@sam_data$inter_rank)
d@sam_data$soil_practice <- gsub("Weeding", "Grassing", d@sam_data$soil_practice)

d@sam_data$terroir <- gsub("Perrier", "Vergèze", d@sam_data$terroir)

Spores number

Code
d@sam_data$nb_spores <- rowMeans(cbind(as.numeric(d@sam_data$spores_rep1), as.numeric(d@sam_data$spores_rep2),as.numeric(d@sam_data$spores_rep3)),na.rm=TRUE)

Practices

Code
d@sam_data$organic <- ifelse(!d@sam_data$practice %in%
  c("Conventional", "Conversion"),
"Organic", "NonOrganic"
)

Paired samples (roots, spores from the same samples)

Code
d@sam_data$paired_name <-
  apply(as.matrix(d@sam_data[, c(5:16)]), 1, paste,
    collapse = "--"
  )
d@sam_data$paired_name[d@sam_data$terroir == "Vergèze"] <-
  gsub("R", "", d@sam_data$ref_mycea[d@sam_data$terroir == "Vergèze"])
Code
tib_sam_data <- as_tibble(d@sam_data) %>%
  mutate(across(starts_with("Myc_freq"), as.numeric)) %>%
  mutate(across(starts_with("Myc_intensity_colonization"), as.numeric)) %>%
  mutate(across(starts_with("Myc_intensity_root"), as.numeric)) %>%
  mutate(across(starts_with("Arbuscul_richness"), as.numeric)) %>%
  mutate(across(starts_with("Arbuscul_abundance"), as.numeric)) %>%
  mutate(across(starts_with("Vesicle_richness"), as.numeric)) %>%
  mutate(across(starts_with("Vesicle_abundance"), as.numeric)) %>%
  mutate(Myc_freq = rowMeans(select(., starts_with("Myc_freq")), na.rm = TRUE)) %>%
  mutate(Myc_intensity_colonization = rowMeans(select(
    ., starts_with("Myc_intensity_colonization")
  ), na.rm = TRUE)) %>%
  mutate(Myc_intensity_root = rowMeans(select(., starts_with(
    "Myc_intensity_root"
  )), na.rm = TRUE)) %>%
  mutate(Arbuscul_richness = rowMeans(select(., starts_with(
    "Arbuscul_richness"
  )), na.rm = TRUE)) %>%
  mutate(Arbuscul_abundance = rowMeans(select(., starts_with(
    "Arbuscul_abundance"
  )), na.rm = TRUE)) %>%
  mutate(Vesicle_richness = rowMeans(
    select(., starts_with(
      "Vesicle_richness"
    )),
    na.rm = TRUE
  )) %>%
  mutate(Vesicle_abundance = rowMeans(
    select(., starts_with(
      "Vesicle_abundance"
    )),
    na.rm = TRUE
  )) %>% tibble::column_to_rownames(var = "X")

sample_data(d) <- tib_sam_data

Exclude Mycelium

Code
d <- clean_pq(subset_samples(d, compartment != "Mycelium"))
d_vs <- clean_pq(subset_samples(tar_read(d_vs), compartment != "Mycelium"))

Summary of final dataset

Code
summary_plot_pq(d) +
  labs(title = "Summary of final dataset") +
  theme(plot.title = element_text(hjust = 0.5, size = 20, color = "#1e2b4c"))

Summary of bioinformatics pipeline

Global summary

Code
kable(tar_read(track_sequences_samples_clusters))
nb_sequences nb_clusters nb_samples
Paired sequences 8283716 28637 194
Paired sequences without chimera 7777304 9791 194
Paired sequences without chimera and longer than 200bp 7736170 7522 194
ASV denoising 7533009 7522 186
OTU after vsearch reclustering at 97% 7533009 1081 186
OTU after blast filter without reclustering 3819912 3793 182
OTU after blast filter with reclustering 3822304 219 182

SUBSETING AND FILTERING DATASET

Code
d_muco <- clean_pq(subset_taxa(d, Class_PR2 == "Mucoromycota"))
d_blast <- clean_pq(subset_taxa(tar_read(d_blast), Class_PR2 == "Mucoromycota"))
Code
for(i in c(19:46)){
  d_muco@sam_data[,i] <- as.numeric(unlist(d_muco@sam_data[,i]))
}
Code
d_roots <- clean_pq(subset_samples(d_muco, compartment == "Roots"))
d_spores <- clean_pq(subset_samples(d_muco, compartment == "Spores"))
Code
d_rs_merged_samples <- clean_pq(merge_samples2(
  d_muco,
  "paired_name",
  default_fun = 
    function(x){MiscMetabar::diff_fct_diff_class(x, na.rm = T)}
))
Code
d_rs <- clean_pq(subset_samples_pq(
    d_rs_merged_samples,
    sample_sums(d_rs_merged_samples) > 2000
  ))

TAXONOMICAL ANALYSES

Code
# Figure S8
multitax_bar_pq(d_vs, "Supergroup_PR2",
  "Subdivision_PR2", "Class_PR2",
  nb_seq = TRUE
) +
  ggtitle("Number of sequences (log10) including non-Mucoromycota")

Number of sequences

Code
# Figure S6
dir.create("output/krona")
Warning in dir.create("output/krona"): 'output/krona' existe déjà
Code
krona(
  d_vs,
  ranks = c(11:19),
  "output/krona/OTU_post_clustered_allEuk",
  name = "OTU_post_clustered_allEuk"
)
krona(
  d,
  ranks = c(11:19),
  "output/krona/OTU_filtOnMaarjam_AM",
  name = "OTU_filtOnMaarjam_AM"
)
krona(
  d_muco,
  ranks = c(11:19),
  "output/krona/OTU_filtOnMaarjamAndPr2_AM",
  name = "OTU_filtOnMaarjamAndPr2_AM"
)
merge_krona(
  c(
    "output/krona/OTU_post_clustered_allEuk",
    "output/krona/OTU_filtOnMaarjam_AM",
    "output/krona/OTU_filtOnMaarjamAndPr2_AM"
  ),
  "output/krona/Euk_to_AM_filtering_nb_seq.html"
)

Number of OTU

Code
# Figure S7
krona(
  d_vs,
  ranks = c(11:19),
  "output/krona/OTU_post_clustered_allEuk_Nbotu",
      name = "OTU_post_clustered_allEuk",
  nb_seq = F
)
krona(
  d,
  ranks = c(11:19),
  "output/krona/OTU_filtOnMaarjam_AM_Nbotu",
    name = "OTU_filtOnMaarjam_AM",
  nb_seq = F
)
krona(
  d_muco,
  ranks = c(11:19),
  "output/krona/OTU_filtOnMaarjamAndPr2_AM_Nbotu",
  name = "OTU_filtOnMaarjamAndPr2_AM",
  nb_seq = F
)
merge_krona(
  c(
    "output/krona/OTU_post_clustered_allEuk_Nbotu",
    "output/krona/OTU_filtOnMaarjam_AM_Nbotu",
    "output/krona/OTU_filtOnMaarjamAndPr2_AM_Nbotu"
  ),
  "output/krona/Euk_to_AM_filtering_Nbotu.html"
)

NUMBER OF SPORES ANALYSES

Code
psm_otu <- psmelt(as_binary_otu_table(d_spores)) |>
  filter(Abundance > 0) |>
  group_by(Sample) |>
  summarise(
    "Abundance" = sum(Abundance),
    "region" = region[1],
    "nb_spores" = nb_spores[1]
  )
Code
psm_res <- psmelt_samples_pq(d_spores) %>% 
  mutate(nb_spores_log10 = log10(nb_spores))

Spores count <-> Alpha diversity

Code
# Fig S2
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_spores, rngseed = 221013)) %>% mutate(nb_spores_log10 = log10(nb_spores))

ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Abundance, type = "non-parametric") +
  ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Hill_0, type = "non-parametric") +
  ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Hill_1, type = "non-parametric") +
  ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Hill_2, type = "non-parametric")

Spores count <-> Practice

Code
# Figure S5
ggstatsplot::ggbetweenstats(psm_res, practice, nb_spores_log10, type = "non-parametric") 

Spores count <-> Terroir

Code
# Figure S4
ggstatsplot::ggbetweenstats(
  mutate(psm_res,
    terroir = reorder(psm_res$terroir, psm_res$nb_spores)
  ),
  terroir,
  nb_spores_log10,
  type = "non-parametric",
  centrality.plotting = F,
  package = "ggthemes",
  palette = "stata_economist"
) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

Spores count <-> Mycorrhization rate

Code
ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Myc_freq, type = "non-parametric") +
  ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Myc_intensity_colonization, type = "non-parametric") +
  ggstatsplot::ggscatterstats(psm_res_rarefied, nb_spores_log10, Arbuscul_abundance, type = "non-parametric")

MYCORRHIZATION MEASURE

Mycorrhization rate <-> Alpha diversity

F = Myc_freq

Code
psm_res <- psmelt_samples_pq(d_rs)
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_rs, rngseed = 221013))

# Figure S3a
(ggstatsplot::ggscatterstats(psm_res, Myc_freq, Hill_0, type = "non-parametrique") +
  ggstatsplot::ggscatterstats(psm_res, Myc_freq, Hill_1, type = "non-parametrique") +
  ggstatsplot::ggscatterstats(psm_res, Myc_freq, Hill_2, type = "non-parametrique")) /
  (ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_freq, Hill_0, type = "non-parametrique") +
    ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_freq, Hill_1, type = "non-parametrique") +
    ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_freq, Hill_2, type = "non-parametrique"))

M = Myc_intensity_colonization

Code
# Figure S3b
(ggstatsplot::ggscatterstats(psm_res, Myc_intensity_colonization, Hill_0, type = "non-parametrique") +
  ggstatsplot::ggscatterstats(psm_res, Myc_intensity_colonization, Hill_1, type = "non-parametrique") +
  ggstatsplot::ggscatterstats(psm_res, Myc_intensity_colonization, Hill_2, type = "non-parametrique")) /
  (ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_intensity_colonization, Hill_0, type = "non-parametrique") +
    ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_intensity_colonization, Hill_1, type = "non-parametrique") +
    ggstatsplot::ggscatterstats(psm_res_rarefied, Myc_intensity_colonization, Hill_2, type = "non-parametrique"))

A = Arbuscul_abundance

Code
# Figure S3c
(ggstatsplot::ggscatterstats(psm_res, Arbuscul_abundance, Hill_0, type = "non-parametrique") +
  ggstatsplot::ggscatterstats(psm_res, Arbuscul_abundance, Hill_1, type = "non-parametrique") +
  ggstatsplot::ggscatterstats(psm_res, Arbuscul_abundance, Hill_2, type = "non-parametrique")) /
  (ggstatsplot::ggscatterstats(psm_res_rarefied, Arbuscul_abundance, Hill_0, type = "non-parametrique") +
    ggstatsplot::ggscatterstats(psm_res_rarefied, Arbuscul_abundance, Hill_1, type = "non-parametrique") +
    ggstatsplot::ggscatterstats(psm_res_rarefied, Arbuscul_abundance, Hill_2, type = "non-parametrique"))

Mycorrhization rate <-> terroir

Code
ggstatsplot::ggbetweenstats(psm_res_rarefied, terroir, Myc_freq, type = "non-parametrique", centrality.plotting = F) +
  ggstatsplot::ggbetweenstats(psm_res_rarefied, terroir, Myc_intensity_colonization, type = "non-parametrique", centrality.plotting = F) +
  ggstatsplot::ggbetweenstats(psm_res_rarefied, terroir, Arbuscul_abundance, type = "non-parametrique", centrality.plotting = F)

Mycorrhization rate <-> practice

Code
ggstatsplot::ggbetweenstats(psm_res_rarefied, practice, Myc_freq, type = "non-parametrique", centrality.plotting = F) +
  ggstatsplot::ggbetweenstats(psm_res_rarefied, practice, Myc_intensity_colonization, type = "non-parametrique", centrality.plotting = F) +
  ggstatsplot::ggbetweenstats(psm_res_rarefied, practice, Arbuscul_abundance, type = "non-parametrique", centrality.plotting = F)

ECOLOGICAL ANALYSES : ALPHA DIVERSITY

Compartment

Code
psm_res <- psmelt_samples_pq(d_muco)
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_muco, rngseed = 221013))
Code
ggstatsplot::ggbetweenstats(psm_res_rarefied,
  compartment,
  Hill_0,
  type = "non-parametrique"
) +
  ggstatsplot::ggbetweenstats(psm_res_rarefied,
    compartment,
    Hill_1,
    type = "non-parametrique"
  ) +
  ggstatsplot::ggbetweenstats(psm_res_rarefied,
    compartment,
    Hill_2,
    type = "non-parametrique"
  )

Code
ggplot(psm_res_rarefied, aes(y=Hill_0, x=compartment, color=practice))+
  geom_point() + 
  geom_line(aes(group = paired_name)) +
  facet_wrap(~terroir)

Terroir

Code
psm_res <- psmelt_samples_pq(d_rs)
psm_res_rarefied <- psmelt_samples_pq(rarefy_even_depth(d_rs, rngseed = 221013))
Code
# Figure 4
 (ggstatsplot::ggbetweenstats(
    mutate(
      psm_res_rarefied,
      terroir = reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_0)
    ),
    terroir,
    Hill_0,
    type = "non-parametrique",
    centrality.plotting = F,
    package = "ggthemes",
    palette = "stata_economist",
    point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
    0.9, size = 3, stroke = 0, na.rm = TRUE),
    violin.args = list(width = 0, linewidth = 0)
  )+ 
  scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
                             "#e15968","#784b43","#c200ab","#00b3f0","#953726",
                             "#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_0))),levels(as.factor(psm_res_rarefied$terroir)))])  + coord_flip()) /
 (ggstatsplot::ggbetweenstats(
    mutate(
      psm_res_rarefied,
      terroir = reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_1)
    ),
    terroir,
    Hill_1,
    type = "non-parametrique",
    centrality.plotting = F,
    package = "ggthemes",
    palette = "stata_economist",
    point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
    0.9, size = 3, stroke = 0, na.rm = TRUE),
    violin.args = list(width = 0, linewidth = 0)
  )+ 
  scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
                             "#e15968","#784b43","#c200ab","#00b3f0","#953726",
                             "#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_1))),levels(as.factor(psm_res_rarefied$terroir)))])  + coord_flip()) /
 (ggstatsplot::ggbetweenstats(
    mutate(
      psm_res_rarefied,
      terroir = reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_2)
    ),
    terroir,
    Hill_2,
    type = "non-parametrique",
    centrality.plotting = F,
    package = "ggthemes",
    palette = "stata_economist",
    point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
    0.9, size = 3, stroke = 0, na.rm = TRUE),
    violin.args = list(width = 0, linewidth = 0)
  )+ 
  scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
                             "#e15968","#784b43","#c200ab","#00b3f0","#953726",
                             "#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res_rarefied$terroir, psm_res_rarefied$Hill_2))),levels(as.factor(psm_res_rarefied$terroir)))])  + coord_flip()) + plot_layout(axes = "collect")

Code
# Figure S10
 (ggstatsplot::ggbetweenstats(
    mutate(
      psm_res,
      terroir = reorder(psm_res$terroir, psm_res$Hill_0)
    ),
    terroir,
    Hill_0,
    type = "non-parametrique",
    centrality.plotting = F,
    package = "ggthemes",
    palette = "stata_economist",
    point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
    0.9, size = 3, stroke = 0, na.rm = TRUE),
    violin.args = list(width = 0, linewidth = 0)
  )+ 
  scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
                             "#e15968","#784b43","#c200ab","#00b3f0","#953726",
                             "#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res$terroir, psm_res$Hill_0))),levels(as.factor(psm_res$terroir)))])  + coord_flip()) /
 (ggstatsplot::ggbetweenstats(
    mutate(
      psm_res,
      terroir = reorder(psm_res$terroir, psm_res$Hill_1)
    ),
    terroir,
    Hill_1,
    type = "non-parametrique",
    centrality.plotting = F,
    package = "ggthemes",
    palette = "stata_economist",
    point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
    0.9, size = 3, stroke = 0, na.rm = TRUE),
    violin.args = list(width = 0, linewidth = 0)
  )+ 
  scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
                             "#e15968","#784b43","#c200ab","#00b3f0","#953726",
                             "#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res$terroir, psm_res$Hill_1))),levels(as.factor(psm_res$terroir)))])  + coord_flip()) /
 (ggstatsplot::ggbetweenstats(
    mutate(
      psm_res,
      terroir = reorder(psm_res$terroir, psm_res$Hill_2)
    ),
    terroir,
    Hill_2,
    type = "non-parametrique",
    centrality.plotting = F,
    package = "ggthemes",
    palette = "stata_economist",
    point.args = list(position = ggplot2::position_jitterdodge(dodge.width = 0.6), alpha =
    0.9, size = 3, stroke = 0, na.rm = TRUE),
    violin.args = list(width = 0, linewidth = 0)
  )+ 
  scale_color_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
                             "#e15968","#784b43","#c200ab","#00b3f0","#953726",
                             "#e09199","#379d30","#ed5019","#ff63b0") [match( levels(as.factor(reorder(psm_res$terroir, psm_res$Hill_2))),levels(as.factor(psm_res$terroir)))])  + coord_flip()) + plot_layout(axes = "collect")

Code
# Fig_7a
p_accu_balanced_mod_terroir <- 
  MiscMetabar::accu_plot_balanced_modality(d_rs,
                                           "terroir", 
                                           nperm = 999, 
                                           step = 300,
                                           quantile_prob=0.9,
                                           sample.size = 10000,
                                           progress_bar = FALSE)
Warning in max(dff$xlab, na.rm = TRUE): aucun argument pour max ; -Inf est
renvoyé
Code
p_accu_balanced_mod_terroir +
  xlab("Number of sequences") + 
  ylab("Number of OTUs") + 
  ggrepel::geom_label_repel(data=summarise(group_by(p_accu_balanced_mod_terroir$data,factor) ,y=max(X1,na.rm = T)),aes(label=factor,y=y,x=max(p_accu_balanced_mod_terroir$data$x,na.rm = T)), max.overlaps = 1900, force = 20) + xlim(c(1,25000)) + theme(legend.position = "none" )
Warning: Removed 14 rows containing missing values or values outside the scale range
(`geom_line()`).

Practice

Code
# Figure 5
ggstatsplot::ggbetweenstats(psm_res_rarefied,
  practice,
  Hill_0,
  centrality.plotting = F,
  type = "non-parametrique"
) + ylab("Hill 0 (richness)") +
  ggstatsplot::ggbetweenstats(psm_res_rarefied,
    practice,
    Hill_1,
    centrality.plotting = F,
    type = "non-parametrique"
  ) + ylab("Hill 1 (~Shannon)") +
  ggstatsplot::ggbetweenstats(psm_res_rarefied,
    practice,
    Hill_2,
    centrality.plotting = F,
    type = "non-parametrique"
  ) + ylab("Hill 2 (~Simpson)") +  plot_layout(axes = "collect")

Code
# Figure S11
ggstatsplot::ggbetweenstats(psm_res,
  practice,
  Hill_0,
  centrality.plotting = F,
  type = "non-parametrique"
) +
  ggstatsplot::ggbetweenstats(psm_res,
    practice,
    Hill_1,
    centrality.plotting = F,
    type = "non-parametrique"
  ) +
  ggstatsplot::ggbetweenstats(psm_res,
    practice,
    Hill_2,
    centrality.plotting = F,
    type = "non-parametrique"
  )

Code
# Fig_7b
p_accu_balanced_mod <- 
  MiscMetabar::accu_plot_balanced_modality(d_rs,
                                           "practice", 
                                           nperm = 999, 
                                           step = 300,
                                           quantile_prob=0.90,
                                           progress_bar = FALSE)
Warning in max(dff$xlab, na.rm = TRUE): aucun argument pour max ; -Inf est
renvoyé
Code
p_accu_balanced_mod +
  xlab("Number of sequences") +
  ylab("Number of OTUs") + 
  ggrepel::geom_label_repel(data=summarise(group_by(p_accu_balanced_mod$data,factor) ,y=max(X1,na.rm = T)),aes(label=factor,y=y,x=max(p_accu_balanced_mod$data$x,na.rm = T)), max.overlaps = 1900, force = 20) + xlim(c(1,18000)) + theme(legend.position = "none" )
Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_line()`).

ZOOM ON LOW-DIVERSITY SAMPLES

Code
otu_hill <- vegan::renyi(d_rs@otu_table, scale = c(0, 1, 2), hill = TRUE)

low_div_samples <-
  rownames(otu_hill[otu_hill$`0` < 10 |
    otu_hill$`1` < 2 |
    otu_hill$`2` < 2, ])

low_div_samples
[1] "PE13"                                                                                                                                           
[2] "2020--Bordelais--Langoiran--Langoiran--Conventional chimique--Conventional---0.356503--44.713204--Conventional--Herbicide--Herbicide--Herbicide"
Code
sum(otu_hill$`0` < 10)
[1] 2
Code
sum(otu_hill$`1` < 2)
[1] 0
Code
sum(otu_hill$`2` < 2)
[1] 0
Code
length(low_div_samples)
[1] 2
Code
low_div_samples_table <- tibble(cbind(as_tibble(d_rs@sam_data[low_div_samples, c("terroir", "practice", "rank", "inter_rank", "compartment")]), apply(otu_hill[low_div_samples, ], 2, round, digits = 2), sample_sums(d_rs)[low_div_samples]))
Warning in class(x) <- tibble_class: La class(x) est une chaîne multiple
("tbl_df", "tbl", ...) ; le résultat ne sera plus un objet S4
Code
low_div_samples_table$compartment <- ifelse(is.na(low_div_samples_table$compartment), "Spores + Roots", "Spores")

colnames(low_div_samples_table) <- c("Terroir", "Global practice", "Rank practice", "Inter rank practice", "Compartment", "Hill 0 (richness)", "Hill 1 (~Shannon)", "Hill 2 (~Simpson)", "nb_seq")

low_div_samples_table <- arrange(low_div_samples_table, as.numeric(`Hill 0 (richness)`))

low_div_samples_table <-
  rbind(
    low_div_samples_table,
    c(
      "",
      "",
      "",
      "MEAN",
      "(low -diversity samples)",
      round(mean(low_div_samples_table$`Hill 0 (richness)`), 2),
      round(mean(low_div_samples_table$`Hill 1 (~Shannon)`), 2),
      round(mean(low_div_samples_table$`Hill 2 (~Simpson)`), 2),
      round(mean(low_div_samples_table$nb_seq), 2)
    ),
    c(
      "",
      "",
      "",
      "MEAN",
      "(all samples)",
      round(mean(otu_hill$`0`), 2),
      round(mean(otu_hill$`1`), 2),
      round(mean(otu_hill$`2`), 2),
      round(mean(sample_sums(d_rs)), 2)
    ),
     c(
      "",
      "",
      "",
      "MAX",
      "(all samples)",
      round(max(otu_hill$`0`), 2),
      round(max(otu_hill$`1`), 2),
      round(max(otu_hill$`2`), 2),
      round(max(sample_sums(d_rs)), 2)
    )
  )

# Table 3
kbl(low_div_samples_table) |>
  kable_classic(full_width = F, html_font = "Cambria")
Terroir Global practice Rank practice Inter rank practice Compartment Hill 0 (richness) Hill 1 (~Shannon) Hill 2 (~Simpson) nb_seq
Langoiran Conventional Herbicide Herbicide Spores + Roots 6 3.05 2.11 12643
Vergèze Organic Scratching Ploughing Spores 7 4.86 3.95 5938
MEAN (low -diversity samples) 6.5 3.96 3.03 9290.5
MEAN (all samples) 101.28 20.09 10.46 48265.19
MAX (all samples) 432 64.2 24.15 101104

ECOLOGICAL ANALYSES : BETA DIVERSITY

Compartment

Code
# Figure S9
p <- phyloseq::plot_ordination(d_muco,
  vegan::decorana(vegdist(as(otu_table(d_muco), "matrix"),
    method = "robust.aitchison"
  )),
  color = "region",
  shape = "compartment"
) +
  geom_point(size = 3) +
  stat_ellipse(inherit.aes = F, aes(x = DCA1, y = DCA2, linetype = compartment))
Warning in phyloseq::plot_ordination(d_muco,
vegan::decorana(vegdist(as(otu_table(d_muco), : `Ordination species/OTU/taxa
coordinate indices did not match `physeq` index names. Setting corresponding
coordinates to NULL.
Code
p + geom_line(data = p$data, aes(group = paired_name))

Practice and terroir

Code
p <-
  phyloseq::plot_ordination(d_rs,
    ordinate(d_rs,
      method = "NMDS",
      distance = "bray"
    ),
    color = "terroir",
    shape = "practice"
  ) + facet_wrap("region") + scale_shape_manual(values=c(25,23,24))
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2199938 
Run 1 stress 0.2203704 
... Procrustes: rmse 0.03925079  max resid 0.2885775 
Run 2 stress 0.2265037 
Run 3 stress 0.2213033 
Run 4 stress 0.2203289 
... Procrustes: rmse 0.03921461  max resid 0.2889235 
Run 5 stress 0.2201826 
... Procrustes: rmse 0.07042354  max resid 0.5261175 
Run 6 stress 0.242704 
Run 7 stress 0.2304432 
Run 8 stress 0.220005 
... Procrustes: rmse 0.07747605  max resid 0.510776 
Run 9 stress 0.2200585 
... Procrustes: rmse 0.01712057  max resid 0.1278995 
Run 10 stress 0.2209168 
Run 11 stress 0.2241108 
Run 12 stress 0.2203529 
... Procrustes: rmse 0.0391036  max resid 0.2892524 
Run 13 stress 0.2192479 
... New best solution
... Procrustes: rmse 0.07702747  max resid 0.514501 
Run 14 stress 0.2196043 
... Procrustes: rmse 0.03466159  max resid 0.2841357 
Run 15 stress 0.2214039 
Run 16 stress 0.2216717 
Run 17 stress 0.2447455 
Run 18 stress 0.2183676 
... New best solution
... Procrustes: rmse 0.03750409  max resid 0.2833134 
Run 19 stress 0.2217418 
Run 20 stress 0.2201252 
*** Best solution was not repeated -- monoMDS stopping criteria:
     1: no. of iterations >= maxit
    19: stress ratio > sratmax
Code
p_nmds <- p + geom_point(
  data = select(p$data, -c(region)),
  fill = "grey40",
  color = "grey20",
  size = 1,
  alpha = 0.5
) +
  geom_point(size = 3, aes(fill=terroir)) +
  ggtitle("NMDS on bray distance") + 
  scale_fill_manual(values=c("#5660d0","#de8a00","#b59e00","#00773d","#fcb709",
                             "#e15968","#784b43","#c200ab","#00b3f0","#953726",
                             "#e09199","#379d30","#ed5019","#ff63b0")) +
  scale_color_manual(values=rep("grey20",14))
# Figure 9
p_nmds

Code
# Figure 6a
upset_pq(d_rs, "practice", taxa_fill = "Family",
         set_sizes=(
        upset_set_size()
         + geom_text(aes(label=after_stat(count)), hjust=-0.3, color="white", stat='count')
    ))

Code
ggvenn_pq(d_rs, "practice")

Code
# Figure 6b
upset_pq(d_rs, "terroir", taxa_fill = "Family", min_size = 2,  height_ratio = 0.6
 ,set_sizes=(
        upset_set_size()
         + geom_label(aes(label=after_stat(count)), hjust=-0.3, size=2.5,  stat='count')
    ))
Warning: Removed 233 rows containing non-finite outside the scale range
(`stat_count()`).

Global beta-div analysis

Spatial data

Code
library(geodist)
lat_lon <- as_tibble(d_rs@sam_data[, c("lat", "long")])
lat_lon$lat <- as.numeric(lat_lon$lat)
lat_lon$long <- as.numeric(lat_lon$long)
dist_spatial_meter <- as.dist(geodist(lat_lon, measure = "geodesic"),
  upper = FALSE
)
Code
lon_lat_rs <- d_rs@sam_data[, c("long", "lat")]
lon_lat_rs$long <- as.numeric(lon_lat_rs$long)
lon_lat_rs$lat <- as.numeric(lon_lat_rs$lat)
MEM <- dbmem(lon_lat_rs, MEM.autocor = "non-null")
Code
test_MEM <- moran.randtest(MEM, nrepet = 999)
test_MEM$pvalue_adjust <- p.adjust(test_MEM$pvalue, method = "BH")
Code
d_rs@sam_data$MEM_1 <- MEM[, 1]
d_rs@sam_data$MEM_2 <- MEM[, 2]
res_ado_spatial_robAit <- adonis_pq(
  d_rs,
  "MEM_1 + MEM_2 + practice + inter_rank + rank  + terroir",
  correction_for_sample_size = TRUE,
  dist_method = "robust.aitchison"
)

res_ado_spatial_robAit_rarefy <- adonis_pq(
  rarefy_even_depth(d_rs, rngseed = 626),
  "MEM_1 + MEM_2 + practice + inter_rank + rank  + terroir",
  correction_for_sample_size = FALSE,
  dist_method = "robust.aitchison"
)

res_ado_spatial_bray <- adonis_pq(
  d_rs,
  "MEM_1 + MEM_2 + practice + inter_rank + rank  + terroir",
  correction_for_sample_size = TRUE,
  dist_method = "bray"
)


res_ado_spatial_bray_rarefy <- adonis_pq(
  rarefy_even_depth(d_rs, rngseed = 626),
  "MEM_1 + MEM_2 + practice + inter_rank + rank  + terroir",
  correction_for_sample_size = FALSE,
  dist_method = "bray"
)

df <- data.frame(res_ado_spatial_bray)
df$names <-
  factor(
    rownames(df),
    levels = c(
      "sample_size",
      "MEM_1",
      "MEM_2",
      "practice",
      "inter_rank",
      "rank",
      "terroir",
      "Residual",
      "Total"
    )
  )
df <- df %>% filter(!names == "Total")
p1 <-
  ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
  geom_bar(stat = "identity") +
  scale_fill_viridis_c() +
  geom_text(
    label = case_when(
      df$Pr..F. < 0.001 ~ "***",
      df$Pr..F. < 0.01 ~ "**",
      df$Pr..F. < 0.05 ~ "*",
      df$Pr..F. > 0.05 ~ "ns"
    ),
    nudge_y = 0.01
  ) +
  ggtitle("Permanova on Bray")

df <- data.frame(res_ado_spatial_robAit)
df$names <-
  factor(
    rownames(df),
    levels = c(
      "MEM_1",
      "MEM_2",
      "practice",
      "inter_rank",
      "rank",
      "terroir",
      "Residual",
      "Total"
    )
  )
df <- df %>% filter(!names == "Total")
p2 <-
  ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
  geom_bar(stat = "identity") +
  scale_fill_viridis_c() +
  geom_text(
    label = case_when(
      df$Pr..F. < 0.001 ~ "***",
      df$Pr..F. < 0.01 ~ "**",
      df$Pr..F. < 0.05 ~ "*",
      df$Pr..F. > 0.05 ~ "ns"
    ),
    nudge_y = 0.01
  ) +
  ggtitle("Permanova on robust Aitchison")

df <- data.frame(res_ado_spatial_bray_rarefy)
df$names <-
  factor(
    rownames(df),
    levels = c(
      "MEM_1",
      "MEM_2",
      "practice",
      "inter_rank",
      "rank",
      "terroir",
      "Residual",
      "Total"
    )
  )
df <- df %>% filter(!names == "Total")

p3 <-
  ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
  geom_bar(stat = "identity") +
  scale_fill_viridis_c() +
  geom_text(
    label = case_when(
      df$Pr..F. < 0.001 ~ "***",
      df$Pr..F. < 0.01 ~ "**",
      df$Pr..F. < 0.05 ~ "*",
      df$Pr..F. > 0.05 ~ "ns"
    ),
    nudge_y = 0.01
  ) +
  ggtitle("Permanova on Bray after rarefaction")

df <- data.frame(res_ado_spatial_robAit_rarefy)
df$names <-
  factor(
    rownames(df),
    levels = c(
      "MEM_1",
      "MEM_2",
      "practice",
      "inter_rank",
      "rank",
      "terroir",
      "Residual",
      "Total"
    )
  )
df <- df %>% filter(!names == "Total")
p4 <-
  ggplot(df, aes(y = R2 / Df, x = names, fill = R2 / Df)) +
  geom_bar(stat = "identity") +
  scale_fill_viridis_c() +
  geom_text(
    label = case_when(
      df$Pr..F. < 0.001 ~ "***",
      df$Pr..F. < 0.01 ~ "**",
      df$Pr..F. < 0.05 ~ "*",
      df$Pr..F. > 0.05 ~ "ns"
    ),
    nudge_y = 0.01
  ) +
  ggtitle("Permanova on robust Aitchison after rarefaction")

(p1 + ylim(0, 0.12) + p2 + ylim(0, 0.12)) /
  (p3 + ylim(0, 0.12) + p4 + ylim(0, 0.12))
Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).
Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).
Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).
Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).

Code
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), d_rs@sam_data$terroir))
Analysis of Variance Table

Response: Distances
          Df  Sum Sq  Mean Sq F value Pr(>F)
Groups    13 0.23494 0.018072  1.4539 0.1618
Residuals 61 0.75823 0.012430               
Code
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), d_rs@sam_data$practice))
Analysis of Variance Table

Response: Distances
          Df  Sum Sq   Mean Sq F value Pr(>F)
Groups     2 0.01739 0.0086972  0.8898 0.4152
Residuals 72 0.70378 0.0097748               
Code
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$terroir))
Analysis of Variance Table

Response: Distances
          Df  Sum Sq  Mean Sq F value Pr(>F)
Groups    13 0.23494 0.018072  1.4539 0.1618
Residuals 61 0.75823 0.012430               
Code
anova(vegan::betadisper(phyloseq::distance(d_rs@otu_table, "bray"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$practice))
Analysis of Variance Table

Response: Distances
          Df  Sum Sq   Mean Sq F value Pr(>F)
Groups     2 0.01739 0.0086972  0.8898 0.4152
Residuals 72 0.70378 0.0097748               
Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$terroir))
Analysis of Variance Table

Response: Distances
          Df  Sum Sq Mean Sq F value    Pr(>F)    
Groups    13 2471.66 190.127  13.793 1.288e-13 ***
Residuals 61  840.86  13.785                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$practice))
Analysis of Variance Table

Response: Distances
          Df Sum Sq Mean Sq F value Pr(>F)
Groups     2   66.0  32.985  0.4173 0.6604
Residuals 72 5691.6  79.050               
Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$terroir))
Analysis of Variance Table

Response: Distances
          Df  Sum Sq Mean Sq F value    Pr(>F)    
Groups    13 2471.66 190.127  13.793 1.288e-13 ***
Residuals 61  840.86  13.785                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), rarefy_even_depth(d_rs, rngseed = 626)@sam_data$practice))
Analysis of Variance Table

Response: Distances
          Df Sum Sq Mean Sq F value Pr(>F)
Groups     2   66.0  32.985  0.4173 0.6604
Residuals 72 5691.6  79.050               

Soil properties

Code
soil_prop <- as_tibble(d_rs@sam_data) |>
  select(paired_name, practice, organic, terroir, Coarse_sand:CEC)
soil_prop_num <- soil_prop |>
  select(-all_of(c("practice", "organic", "terroir", "CaO", "Total_sand", "Total_filt"))) |>
  tibble::column_to_rownames("paired_name") |>
  mutate(across(everything(), as.numeric)) |>
  tidyr::drop_na()
Code
pca_test_res <- PCAtest::PCAtest(soil_prop_num, varcorr = T, plot = F, counter = F)

Sampling bootstrap replicates... Please wait

Calculating confidence intervals of empirical statistics... Please wait

Sampling random permutations... Please wait

Comparing empirical statistics with their null distributions... Please wait

========================================================
Test of PCA significance: 17 variables, 44 observations
1000 bootstrap replicates, 1000 random permutations
========================================================

Empirical Psi = 55.8085, Max null Psi = 9.1059, Min null Psi = 4.3334, p-value = 0
Empirical Phi = 0.4530, Max null Phi = 0.1830, Min null Phi = 0.1262, p-value = 0

Empirical eigenvalue #1 = 7.28487, Max null eigenvalue = 3.01024, p-value = 0
Empirical eigenvalue #2 = 3.15651, Max null eigenvalue = 2.34796, p-value = 0
Empirical eigenvalue #3 = 2.59627, Max null eigenvalue = 2.04867, p-value = 0
Empirical eigenvalue #4 = 1.16046, Max null eigenvalue = 1.85095, p-value = 1
Empirical eigenvalue #5 = 0.93253, Max null eigenvalue = 1.63065, p-value = 1
Empirical eigenvalue #6 = 0.74202, Max null eigenvalue = 1.51873, p-value = 1
Empirical eigenvalue #7 = 0.39204, Max null eigenvalue = 1.30991, p-value = 1
Empirical eigenvalue #8 = 0.24218, Max null eigenvalue = 1.18896, p-value = 1
Empirical eigenvalue #9 = 0.19075, Max null eigenvalue = 1.04599, p-value = 1
Empirical eigenvalue #10 = 0.08962, Max null eigenvalue = 0.95824, p-value = 1
Empirical eigenvalue #11 = 0.07595, Max null eigenvalue = 0.83144, p-value = 1
Empirical eigenvalue #12 = 0.05693, Max null eigenvalue = 0.76089, p-value = 1
Empirical eigenvalue #13 = 0.0369, Max null eigenvalue = 0.65485, p-value = 1
Empirical eigenvalue #14 = 0.02507, Max null eigenvalue = 0.56529, p-value = 1
Empirical eigenvalue #15 = 0.01767, Max null eigenvalue = 0.51694, p-value = 1
Empirical eigenvalue #16 = 0.00018, Max null eigenvalue = 0.3948, p-value = 1
Empirical eigenvalue #17 = 5e-05, Max null eigenvalue = 0.30934, p-value = 1

PC 1 is significant and accounts for 42.9% (95%-CI:39.1-49.6) of the total variation
PC 2 is significant and accounts for 18.6% (95%-CI:16.7-24.4) of the total variation
PC 3 is significant and accounts for 15.3% (95%-CI:10.2-17.9) of the total variation

The first 3 PC axes are significant and account for 76.7% of the total variation

Variables 1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 12, 13, 15, 16, and 17 have significant loadings on PC 1
Variables 4, 8, 10, 12, and 14 have significant loadings on PC 2
Variables 5, 6, and 7 have significant loadings on PC 3

Variables 8, 10, 13, 15, and 17 have significant correlations with PC 1
Variables 8, 10, 12, and 14 have significant correlations with PC 2
Variables , and  have significant correlations with PC 3
Code
pval <- c()
for (i in seq_len(length(pca_test_res$`Percentage of variation of empirical PC's`))) {
  obs <- pca_test_res$`Percentage of variation of empirical PC's`[i]
  null_model <- pca_test_res$`Percentage of variation of randomized data`[, i]

  pval[i] <- (sum(obs < null_model) + 1) / (1 + nrow(pca_test_res$`Percentage of variation of randomized data`))
}

pca_test_pval_adj <- p.adjust(pval, method = "BH")
pca_test_pval_adj
 [1] 0.005661006 0.005661006 0.005661006 1.000000000 1.000000000 1.000000000
 [7] 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000
[13] 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000
Code
res_pca <- dudi.pca(soil_prop_num, scannf = F, nf = 4)
Code
p_pca_variable <- fviz_pca_var(res_pca,
  col.var = "cos2",
  gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
  repel = TRUE
)

p_pca_variable_axe3_4 <- fviz_pca_var(res_pca,
  col.var = "cos2", axes = c(3, 4),
  gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
  repel = TRUE
)

res_pca_var <- get_pca_var(res_pca)

# Fig. 2
(p_pca_variable + p_pca_variable_axe3_4) / free((factoextra::fviz_eig(res_pca)) +
  (ggcorrplot::ggcorrplot(res_pca_var$cor)))

Code
practice <- unlist(lapply(strsplit(rownames(res_pca$tab), "--"), function(x) {
  x[5]
}))

terroir <- unlist(lapply(strsplit(rownames(res_pca$tab), "--"), function(x) {
  x[2]
}))

p_pca_ind_terroir <- fviz_pca_ind(res_pca,
  habillage = terroir,
  label = FALSE,
  repel = TRUE,
  mean.point = FALSE,
  addEllipses = TRUE,
  pointsize = 3,
  ellipse.type = "convex"
) +
  scale_shape_manual(values = rep(16, length(unique(terroir))))

p_pca_ind_terroir$data$terroir <- terroir
p_pca_ind_terroir$data$practice <- practice

p_pca_ind_terroir_practice <- p_pca_ind_terroir + geom_point(aes(x = x, y = y),
  shape = case_when(
    practice == "Organic" ~ 2,
    practice == "Conventional" ~ 6,
    practice == "Conversion" ~ 5
  ), size = 3.5,
  color = rgb(0, 0, 0, 0.7)
)


# value
kw1_terroir <- kruskal.test(res_pca$tab[, 1], terroir)
kw2_terroir <- kruskal.test(res_pca$tab[, 2], terroir)
kw3_terroir <- kruskal.test(res_pca$tab[, 3], terroir)

kw1_practice <- kruskal.test(res_pca$tab[, 1], practice)
kw2_practice <- kruskal.test(res_pca$tab[, 2], practice)
kw3_practice <- kruskal.test(res_pca$tab[, 3], practice)

# Figure 3
p_pca_ind_terroir_practice + 
  annotate("text", x = 2, y = 4.5, label = paste("Kruskal-Wallis (axe1 - terroir): pval=", round(as.numeric(kw1_terroir$p.value),4)), size = 3) + 
  annotate("text", x = 2, y = 4.25, label = paste("Kruskal-Wallis (axe2 - terroir): pval=", round(as.numeric(kw2_terroir$p.value),4)), size = 3) + 
  annotate("text", x = 2, y = 4, label = paste("Kruskal-Wallis (axe3 - terroir): pval=", round(as.numeric(kw3_terroir$p.value),4)), size = 3) + 
  annotate("text", x = -2, y = 4.5, label = paste("Kruskal-Wallis (axe1 - practice): pval=", round(as.numeric(kw1_practice$p.value),4)), size = 3) + 
  annotate("text", x = -2, y = 4.25, label = paste("Kruskal-Wallis (axe2 - practice): pval=", round(as.numeric(kw2_practice$p.value),4)), size = 3) + 
  annotate("text", x = -2, y = 4, label = paste("Kruskal-Wallis (axe3 - practice): pval=", round(as.numeric(kw3_practice$p.value),4)), size = 3)
Warning: Removed 21 rows containing missing values or values outside the scale range
(`geom_point()`).

Code
# Figure S1
ggplot(d_rs@sam_data, aes(x = practice, y = as.numeric(Cu), fill = practice)) +
  geom_violin() +
  geom_jitter()
Warning: Removed 23 rows containing non-finite outside the scale range
(`stat_ydensity()`).
Warning: Removed 23 rows containing missing values or values outside the scale range
(`geom_point()`).

Code
kruskal.test(as.numeric(d_rs@sam_data$Cu), d_rs@sam_data$practice)

    Kruskal-Wallis rank sum test

data:  as.numeric(d_rs@sam_data$Cu) and d_rs@sam_data$practice
Kruskal-Wallis chi-squared = 0.72149, df = 2, p-value = 0.6972
Code
d_with_pca <-
  clean_pq(subset_samples_pq(d_rs, d_rs@sam_data$paired_name %in% rownames(soil_prop_num)))

res_pca_ind <- get_pca_ind(res_pca)

if (!identical(match(rownames(res_pca_ind$coord), d_with_pca@sam_data$paired_name), sort(match(rownames(res_pca_ind$coord), d_with_pca@sam_data$paired_name)))) {
  stop("ERROR")
}

d_with_pca <- add_info_to_sam_data(d_with_pca, res_pca_ind$coord)
Code
dim_df <- as_tibble(d_with_pca@sam_data) |>
  select(c(starts_with("Dim"), "nb_seq", "nb_otu"))

cor_results <- correlation::correlation(dim_df)
cor_results %>%
  summary(redundant = TRUE) %>%
  plot()

Code
lon_lat_with_pca <- d_with_pca@sam_data[, c("long", "lat")]
lon_lat_with_pca$long <- as.numeric(lon_lat_with_pca$long)
lon_lat_with_pca$lat <- as.numeric(lon_lat_with_pca$lat)
MEM_with_pca <- dbmem(lon_lat_with_pca, MEM.autocor = "non-null")

d_with_pca@sam_data$MEM_1 <- MEM_with_pca[, 1]
d_with_pca@sam_data$MEM_2 <- MEM_with_pca[, 2]
Code
res_ado_spatial_soil_bray <-
  adonis_pq(
    d_with_pca,
    "MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
    correction_for_sample_size = TRUE
  )
res_ado_spatial_soil_robAit <-
  adonis_pq(
    d_with_pca,
    "MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
    correction_for_sample_size = TRUE,
    dist_method = "robust.aitchison"
  )

res_ado_spatial_soil_bray_rarefy <-
  adonis_pq(
    rarefy_even_depth(d_with_pca, rngseed = 626),
    "MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
    correction_for_sample_size = FALSE
  )
res_ado_spatial_soil_robAit_rarefy <-
  adonis_pq(
    rarefy_even_depth(d_with_pca, rngseed = 626),
    "MEM_1 + MEM_2 + Dim.1 + Dim.2 + Dim.3 + practice + inter_rank + rank + terroir",
    correction_for_sample_size = FALSE,
    dist_method = "robust.aitchison"
  )

anova(betadisper(phyloseq::distance(d_rs, method = "bray"), d_rs@sam_data$practice))
Analysis of Variance Table

Response: Distances
          Df  Sum Sq   Mean Sq F value Pr(>F)
Groups     2 0.01739 0.0086972  0.8898 0.4152
Residuals 72 0.70378 0.0097748               
Code
anova(betadisper(phyloseq::distance(d_rs, method = "bray"), d_rs@sam_data$terroir))
Analysis of Variance Table

Response: Distances
          Df  Sum Sq  Mean Sq F value Pr(>F)
Groups    13 0.23494 0.018072  1.4539 0.1618
Residuals 61 0.75823 0.012430               
Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$practice))
Analysis of Variance Table

Response: Distances
          Df Sum Sq Mean Sq F value Pr(>F)
Groups     2   66.0  32.985  0.4173 0.6604
Residuals 72 5691.6  79.050               
Code
anova(vegan::betadisper(vegdist(d_rs@otu_table, "robust.aitchison"), d_rs@sam_data$terroir))
Analysis of Variance Table

Response: Distances
          Df  Sum Sq Mean Sq F value    Pr(>F)    
Groups    13 2471.66 190.127  13.793 1.288e-13 ***
Residuals 61  840.86  13.785                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# TABLE 2
kbl(res_ado_spatial_bray) |>
  kable_classic(full_width = F, html_font = "Cambria")
Df SumOfSqs R2 F Pr(>F)
sample_size 1 0.7692961 0.0353392 3.462441 0.001
MEM_1 1 1.5042988 0.0691030 6.770535 0.001
MEM_2 1 1.1007448 0.0505650 4.954223 0.001
practice 2 0.8165301 0.0375090 1.837516 0.007
inter_rank 3 0.7332177 0.0336819 1.100020 0.285
rank 2 0.6416712 0.0294765 1.444014 0.060
terroir 13 4.8718242 0.2237972 1.686697 0.001
Residual 51 11.3313406 0.5205283 NA NA
Total 74 21.7689236 1.0000000 NA NA
Code
kbl(res_ado_spatial_soil_bray) |>
  kable_classic(full_width = F, html_font = "Cambria")
Df SumOfSqs R2 F Pr(>F)
sample_size 1 0.7801743 0.0603883 3.5267019 0.001
MEM_1 1 0.4296466 0.0332562 1.9421759 0.012
MEM_2 1 1.4246560 0.1102734 6.4400191 0.001
Dim.1 1 0.6061032 0.0469145 2.7398304 0.003
Dim.2 1 0.4590091 0.0355289 2.0749062 0.007
Dim.3 1 0.2666740 0.0206415 1.2054739 0.215
practice 2 0.5569904 0.0431130 1.2589105 0.139
inter_rank 2 0.3033908 0.0234835 0.6857243 0.914
rank 2 0.5008934 0.0387709 1.1321201 0.285
terroir 11 3.1673809 0.2451665 1.3016216 0.016
Residual 20 4.4243845 0.3424631 NA NA
Total 43 12.9193032 1.0000000 NA NA
Code
# TABLE S2
kbl(res_ado_spatial_bray_rarefy) |>
  kable_classic(full_width = F, html_font = "Cambria")
Df SumOfSqs R2 F Pr(>F)
MEM_1 1 1.5265059 0.0715964 6.9137648 0.001
MEM_2 1 1.3436557 0.0630203 6.0856097 0.001
practice 2 0.9547461 0.0447796 2.1620913 0.003
inter_rank 3 0.6228395 0.0292125 0.9403099 0.549
rank 2 0.6161566 0.0288991 1.3953309 0.082
terroir 13 4.7758878 0.2239993 1.6638987 0.001
Residual 52 11.4811988 0.5384928 NA NA
Total 74 21.3209903 1.0000000 NA NA
Code
# TABLE S3
kbl(res_ado_spatial_soil_bray_rarefy) |>
  kable_classic(full_width = F, html_font = "Cambria")
Df SumOfSqs R2 F Pr(>F)
MEM_1 1 0.6248397 0.0496723 2.877821 0.001
MEM_2 1 1.4438602 0.1147813 6.649979 0.001
Dim.1 1 0.6499372 0.0516675 2.993412 0.003
Dim.2 1 0.5208441 0.0414051 2.398849 0.005
Dim.3 1 0.3246755 0.0258104 1.495356 0.088
practice 2 0.5707724 0.0453742 1.314402 0.140
inter_rank 2 0.2656581 0.0211188 0.611770 0.966
rank 2 0.5184269 0.0412129 1.193858 0.194
terroir 11 3.1006425 0.2464891 1.298238 0.030
Residual 21 4.5595725 0.3624684 NA NA
Total 43 12.5792291 1.0000000 NA NA
Code
# TABLE S4
kbl(res_ado_spatial_robAit_rarefy) |>
  kable_classic(full_width = F, html_font = "Cambria")
Df SumOfSqs R2 F Pr(>F)
MEM_1 1 305.5208 0.0226202 1.9770369 0.001
MEM_2 1 277.9747 0.0205808 1.7987853 0.005
practice 2 428.4877 0.0317245 1.3863802 0.002
inter_rank 3 417.9781 0.0309464 0.9015842 0.804
rank 2 433.7451 0.0321138 1.4033908 0.010
terroir 13 3607.0130 0.2670571 1.7954709 0.001
Residual 52 8035.8037 0.5949572 NA NA
Total 74 13506.5231 1.0000000 NA NA
Code
# TABLE S5
kbl(res_ado_spatial_soil_robAit_rarefy) |>
  kable_classic(full_width = F, html_font = "Cambria")
Df SumOfSqs R2 F Pr(>F)
MEM_1 1 433.9694 0.0317449 1.4854760 0.001
MEM_2 1 422.1356 0.0308793 1.4449687 0.002
Dim.1 1 582.6894 0.0426238 1.9945441 0.001
Dim.2 1 392.5284 0.0287135 1.3436235 0.018
Dim.3 1 402.5363 0.0294456 1.3778805 0.021
practice 2 604.4120 0.0442128 1.0344502 0.330
inter_rank 2 444.2780 0.0324990 0.7603811 0.954
rank 2 591.4238 0.0432627 1.0122210 0.424
terroir 11 3661.5622 0.2678439 1.1394106 0.022
Residual 21 6134.9747 0.4487744 NA NA
Total 43 13670.5097 1.0000000 NA NA

Variance partitionning

Code
# Figure 8a
res_varpart_rarefy <- var_par_rarperm_pq(
  physeq = d_with_pca,
  list_component = list(
    "Spatial" = c("MEM_1", "MEM_2"),
    "Soil" = c("Dim.1", "Dim.2", "Dim.3"),
    "Practice" = c("practice", "inter_rank", "rank"),
    "Terroir" = c("terroir")
  ),
  nperm = 99,
  dbrda_computation = TRUE,
  progress_bar = FALSE
)
plot_var_part_pq(res_varpart_rarefy, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)

Code
# Figure 8b
res_varpart_rarefy_wo_soil <- var_par_rarperm_pq(
  physeq = d_rs,
  list_component = list(
    "Spatial" = c("MEM_1", "MEM_2"),
    "Practice" = c("practice", "inter_rank", "rank"),
    "Terroir" = c("terroir")
  ),
  nperm = 99,
  dbrda_computation = TRUE,
  progress_bar = FALSE
)

plot_var_part_pq(res_varpart_rarefy_wo_soil, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)

Code
res_varpart_robAit_rarefy <- var_par_rarperm_pq(
  physeq = d_with_pca,
  list_component = list(
    "Spatial" = c("MEM_1", "MEM_2"),
    "Soil" = c("Dim.1", "Dim.2", "Dim.3"),
    "Practice" = c("practice", "inter_rank", "rank"),
    "Terroir" = c("terroir")
  ),
  dist_method = "robust.aitchison",
  nperm = 99,
  progress_bar = FALSE
)
plot_var_part_pq(res_varpart_robAit_rarefy, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)

Code
# Figure S12b
res_varpart_robAit_rarefy_wo_soil <- var_par_rarperm_pq(
  physeq = d_rs,
  list_component = list(
    "Spatial" = c("MEM_1", "MEM_2"),
    "Practice" = c("practice", "inter_rank", "rank"),
    "Terroir" = c("terroir")
  ),
    dist_method = "robust.aitchison",
  nperm = 99,
  dbrda_computation = TRUE,
  progress_bar = FALSE
)

plot_var_part_pq(res_varpart_robAit_rarefy_wo_soil, show_quantiles = FALSE, filter_quantile_zero = TRUE, show_dbrda_signif = F,digits = 2)

Differential abundance analysis

Indicspecies

Terroir

Code
# Figure 10 
res_mpt_terroir <-
  multipatt(as.matrix(d_rs@otu_table),
    d_rs@sam_data$terroir,
    control = how(nperm = 9999),
    max.order = 3
  )

res_mpt_terroir_df <- res_mpt_terroir$sign
res_mpt_terroir_df$p.adj <- p.adjust(res_mpt_terroir_df$p.value, method = "BH")
res_mpt_terroir_df$OTU_names <- rownames(res_mpt_terroir_df)
res_mpt_terroir_df_signif <-
  res_mpt_terroir_df %>%
  filter(p.adj < 0.05) %>%
  tidyr::pivot_longer(cols = starts_with("s."))

tax_tab <- as.data.frame(d_rs@tax_table)
tax_tab$otu <- rownames(tax_tab)

res_mpt_terroir_df_signif_taxo <- left_join(res_mpt_terroir_df_signif,tax_tab, by = join_by("OTU_names" == "otu"))

ggplot(
  res_mpt_terroir_df_signif_taxo,
  aes(
    x = OTU_names,
    y = name,
    alpha = value,
    color = Genus
  )
) +
  geom_point(size = 6) +
  theme(axis.text.x = element_text(
    angle = 90,
    vjust = 0.5,
    hjust = 1
  ))+
  scale_alpha(range = c(0, 1))

Code
res_mpt_terroir_rg <-
  multipatt(as.matrix(rarefy_even_depth(d_rs, rngseed = 626)@otu_table),
    d_rs@sam_data$terroir,
    control = how(nperm = 9999),
    max.order = 3,
    func = "r.g"
  )

res_mpt_terroir_rg_df <- res_mpt_terroir_rg$sign
res_mpt_terroir_rg_df$p.adj <- p.adjust(res_mpt_terroir_rg_df$p.value, method = "BH")
res_mpt_terroir_rg_df$OTU_names <- rownames(res_mpt_terroir_rg_df)
res_mpt_terroir_rg_df_signif <-
  res_mpt_terroir_rg_df %>%
  filter(p.adj < 0.05) %>%
  tidyr::pivot_longer(cols = starts_with("s."))

ggplot(
  res_mpt_terroir_rg_df_signif,
  aes(
    x = OTU_names,
    y = name,
    alpha = value,
    color = stat
  )
) +
  geom_point(size = 4) +
  theme(axis.text.x = element_text(
    angle = 90,
    vjust = 0.5,
    hjust = 1
  ))

practice

Code
res_mpt_practice <-
  multipatt(as.matrix(d_rs@otu_table),
    d_rs@sam_data$practice,
    control = how(nperm = 9999),
    max.order = 3
  )

res_mpt_practice_df <- res_mpt_practice$sign
res_mpt_practice_df$p.adj <- p.adjust(res_mpt_practice_df$p.value, method = "BH")
res_mpt_practice_df$OTU_names <- rownames(res_mpt_practice_df)
res_mpt_practice_df_signif <-
  res_mpt_practice_df %>%
  filter(p.adj < 0.05) %>%
  tidyr::pivot_longer(cols = starts_with("s."))
res_mpt_practice_df_signif
# A tibble: 0 × 7
# ℹ 7 variables: index <int>, stat <dbl>, p.value <dbl>, p.adj <dbl>,
#   OTU_names <chr>, name <chr>, value <dbl>
Code
res_mpt_practice_rg <-
  multipatt(as.matrix(d_rs@otu_table),
    d_rs@sam_data$practice,
    control = how(nperm = 9999),
    max.order = 3,
    func = "r.g"
  )

res_mpt_practice_rg_df <- res_mpt_practice_rg$sign
res_mpt_practice_rg_df$p.adj <- p.adjust(res_mpt_practice_rg_df$p.value, method = "BH")
res_mpt_practice_rg_df$OTU_names <- rownames(res_mpt_practice_rg_df)
res_mpt_practice_rg_df_signif <-
  res_mpt_practice_rg_df %>%
  filter(p.adj < 0.05) %>%
  tidyr::pivot_longer(cols = starts_with("s."))

res_mpt_practice_rg_df_signif
# A tibble: 0 × 7
# ℹ 7 variables: index <int>, stat <dbl>, p.value <dbl>, p.adj <dbl>,
#   OTU_names <chr>, name <chr>, value <dbl>

Session information

Code
sessionInfo()
R version 4.3.3 (2024-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 11 (bullseye)

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.13.so;  LAPACK version 3.9.0

locale:
 [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    
 [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
 [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/Paris
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] geodist_0.0.8               ComplexUpset_1.3.3         
 [3] microbiomeMarker_1.8.0      formattable_0.2.1          
 [5] kableExtra_1.4.0            adegraphics_1.0-21         
 [7] ade4_1.7-22                 adespatial_0.3-23          
 [9] sf_1.0-16                   indicspecies_1.7.14        
[11] DESeq2_1.42.1               SummarizedExperiment_1.32.0
[13] Biobase_2.62.0              MatrixGenerics_1.14.0      
[15] matrixStats_1.3.0           GenomicRanges_1.54.1       
[17] GenomeInfoDb_1.38.8         IRanges_2.36.0             
[19] S4Vectors_0.40.2            BiocGenerics_0.48.1        
[21] factoextra_1.0.7            FactoMineR_2.11            
[23] vegan_2.6-6.1               lattice_0.22-6             
[25] permute_0.9-7               patchwork_1.2.0            
[27] knitr_1.46                  DT_0.33                    
[29] metacoder_0.3.7             targets_1.7.0              
[31] MiscMetabar_0.9.2           purrr_1.0.2                
[33] dplyr_1.1.4                 dada2_1.30.0               
[35] Rcpp_1.0.12                 ggplot2_3.5.1              
[37] phyloseq_1.46.0            

loaded via a namespace (and not attached):
  [1] igraph_2.0.3                    mia_1.10.0                     
  [3] Formula_1.2-5                   scater_1.30.1                  
  [5] rematch2_2.1.2                  zlibbioc_1.48.2                
  [7] tidyselect_1.2.1                bit_4.0.5                      
  [9] doParallel_1.0.17               BWStest_0.2.3                  
 [11] clue_0.3-65                     rjson_0.2.21                   
 [13] blob_1.2.4                      stringr_1.5.1                  
 [15] rngtools_1.5.2                  S4Arrays_1.2.1                 
 [17] parallel_4.3.3                  RNeXML_2.4.11                  
 [19] correlation_0.8.4               png_0.1-8                      
 [21] cli_3.6.3                       ggplotify_0.1.2                
 [23] CVXR_1.0-12                     bayestestR_0.13.2              
 [25] multtest_2.58.0                 MultiAssayExperiment_1.28.0    
 [27] bluster_1.12.0                  BiocNeighbors_1.20.2           
 [29] TreeSummarizedExperiment_2.10.0 adegenet_2.1.10                
 [31] mime_0.12                       evaluate_0.23                  
 [33] tidytree_0.4.6                  ComplexHeatmap_2.18.0          
 [35] ggh4x_0.2.8                     stringi_1.8.4                  
 [37] backports_1.4.1                 PMCMRplus_1.9.10               
 [39] lmerTest_3.1-3                  gsl_2.1-8                      
 [41] XML_3.99-0.16.1                 Exact_3.2                      
 [43] ggVennDiagram_1.5.2             httpuv_1.6.15                  
 [45] Wrench_1.20.0                   paletteer_1.6.0                
 [47] magrittr_2.0.3                  leaps_3.1                      
 [49] splines_4.3.3                   jpeg_0.1-10                    
 [51] doRNG_1.8.6                     wk_0.9.1                       
 [53] rootSolve_1.8.2.4               ggbeeswarm_0.7.2               
 [55] DBI_1.2.2                       withr_3.0.0                    
 [57] class_7.3-22                    systemfonts_1.0.6              
 [59] htmlwidgets_1.6.4               fs_1.6.4                       
 [61] SuppDists_1.1-9.7               SingleCellExperiment_1.24.0    
 [63] ggrepel_0.9.5                   labeling_0.4.3                 
 [65] SparseArray_1.2.4               cellranger_1.1.0               
 [67] flashClust_1.01-2               rncl_0.8.7                     
 [69] see_0.8.4                       lmom_3.0                       
 [71] effectsize_0.8.8                zoo_1.8-12                     
 [73] XVector_0.42.0                  decontam_1.22.0                
 [75] secretbase_0.5.0                foreach_1.5.2                  
 [77] fansi_1.0.6                     caTools_1.18.2                 
 [79] grid_4.3.3                      data.table_1.15.4              
 [81] ggtree_3.10.1                   rhdf5_2.46.1                   
 [83] irlba_2.3.5.1                   gridGraphics_0.5-1             
 [85] DescTools_0.99.54               base64url_1.4                  
 [87] lazyeval_0.2.2                  yaml_2.3.8                     
 [89] survival_3.6-4                  crayon_1.5.3                   
 [91] RColorBrewer_1.1-3              tidyr_1.3.1                    
 [93] later_1.3.2                     codetools_0.2-19               
 [95] base64enc_0.1-3                 GlobalOptions_0.1.2            
 [97] shape_1.4.6.1                   limma_3.58.1                   
 [99] estimability_1.5.1              Rsamtools_2.18.0               
[101] ggside_0.3.1                    foreign_0.8-86                 
[103] pkgconfig_2.0.3                 ggpubr_0.6.0                   
[105] xml2_1.3.6                      GenomicAlignments_1.38.2       
[107] aplot_0.2.2                     scatterplot3d_0.3-44           
[109] ape_5.8                         viridisLite_0.4.2              
[111] xtable_1.8-4                    interp_1.1-6                   
[113] car_3.1-2                       highr_0.10                     
[115] hwriter_1.3.2.1                 plyr_1.8.9                     
[117] httr_1.4.7                      rbibutils_2.2.16               
[119] tools_4.3.3                     broom_1.0.5                    
[121] beeswarm_0.4.0                  htmlTable_2.4.2                
[123] checkmate_2.3.1                 nlme_3.1-164                   
[125] PCAtest_0.0.1                   lme4_1.1-35.3                  
[127] digest_0.6.36                   numDeriv_2016.8-1.1            
[129] adephylo_1.1-16                 Matrix_1.6-5                   
[131] farver_2.1.2                    reshape2_1.4.4                 
[133] yulab.utils_0.1.4               viridis_0.6.5                  
[135] DirichletMultinomial_1.44.0     rpart_4.1.23                   
[137] glue_1.7.0                      cachem_1.0.8                   
[139] Hmisc_5.1-2                     generics_0.1.3                 
[141] Biostrings_2.70.3               classInt_0.4-10                
[143] mvtnorm_1.2-4                   spdep_1.3-3                    
[145] metagenomeSeq_1.43.0            statmod_1.5.0                  
[147] biomformat_1.30.0               ScaledMatrix_1.10.0            
[149] carData_3.0-5                   minqa_1.2.6                    
[151] ANCOMBC_2.4.0                   glmnet_4.1-8                   
[153] utf8_1.2.4                      seqinr_4.2-36                  
[155] gtools_3.9.5                    readxl_1.4.3                   
[157] datawizard_0.10.0               ggsignif_0.6.4                 
[159] gridExtra_2.3                   shiny_1.8.1.1                  
[161] GenomeInfoDbData_1.2.11         energy_1.7-11                  
[163] rhdf5filters_1.14.1             parameters_0.21.6              
[165] RCurl_1.98-1.14                 memoise_2.0.1                  
[167] rmarkdown_2.26                  scales_1.3.0                   
[169] gld_2.6.6                       svglite_2.1.3                  
[171] phylobase_0.8.12                rstudioapi_0.16.0              
[173] cluster_2.1.6                   hms_1.1.3                      
[175] ggcorrplot_0.1.4.1              munsell_0.5.1                  
[177] colorspace_2.1-0                rlang_1.1.4                    
[179] s2_1.1.6                        DelayedMatrixStats_1.24.0      
[181] sparseMatrixStats_1.14.0        circlize_0.4.16                
[183] scuttle_1.12.0                  mgcv_1.9-1                     
[185] xfun_0.43                       multcompView_0.1-10            
[187] coda_0.19-4.1                   e1071_1.7-14                   
[189] TH.data_1.1-2                   plotROC_2.3.1                  
[191] iterators_1.0.14                emmeans_1.10.1                 
[193] abind_1.4-5                     tibble_3.2.1                   
[195] treeio_1.26.0                   gmp_0.7-4                      
[197] Rhdf5lib_1.24.2                 DECIPHER_2.30.0                
[199] bitops_1.0-7                    Rdpack_2.6                     
[201] ps_1.7.6                        promises_1.3.0                 
[203] RSQLite_2.3.6                   sandwich_3.1-0                 
[205] DelayedArray_0.28.0             proxy_0.4-27                   
[207] Rmpfr_0.9-5                     compiler_4.3.3                 
[209] statsExpressions_1.5.4          prettyunits_1.2.0              
[211] boot_1.3-30                     beachmat_2.18.1                
[213] BiocSingular_1.18.0             units_0.8-5                    
[215] MASS_7.3-60.0.1                 kSamples_1.2-10                
[217] progress_1.2.3                  uuid_1.2-0                     
[219] BiocParallel_1.36.0             insight_0.19.11                
[221] R6_2.5.1                        rstatix_0.7.2                  
[223] fastmap_1.1.1                   multcomp_1.4-25                
[225] vipor_0.4.7                     rsvd_1.0.5                     
[227] ggstatsplot_0.12.3              nnet_7.3-19                    
[229] gtable_0.3.5                    KernSmooth_2.23-22             
[231] latticeExtra_0.6-30             deldir_2.0-4                   
[233] htmltools_0.5.8.1               RcppParallel_5.1.7             
[235] bit64_4.0.5                     lifecycle_1.0.4                
[237] processx_3.8.4                  spData_2.3.0                   
[239] nloptr_2.0.3                    callr_3.7.6                    
[241] vctrs_0.6.5                     zeallot_0.1.0                  
[243] ggfun_0.1.4                     sp_2.1-4                       
[245] pillar_1.9.0                    gplots_3.1.3.1                 
[247] prismatic_1.1.2                 ShortRead_1.60.0               
[249] locfit_1.5-9.9                  jsonlite_1.8.8                 
[251] expm_0.999-9                    GetoptLong_1.0.5